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Abstract. We study the nonequilibrium phase transition in a contact process 
with extended quenched defects by means of Monte-Carlo simulations. We find 
that the spatial disorder correlations dramatically increase the effects of the 
impurities. As a result, the sharp phase transition is completely destroyed by 
smearing. This is caused by effects similar to but stronger than the usual Griffiths 
phenomena, viz., rare strongly coupled spatial regions can undergo the phase 
transition independently from the bulk system. We determine both the stationary 
density in the vicinity of the smeared transition and its time evolution, and we 
compare the simulation results to a recent theory based on extremal statistics. 
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1. Introduction 

Rare regions are an important, if intricate, aspect of systems with impurities and 
defects. In recent years, their influence on phase transitions and critical phenomena 
has reattracted considerable attention. Rare region effects were first studied in the 
context of classical equilibrium phase transitions. Griffiths [1] showed that they lead 
to a singular free energy in an entire parameter region in the vicinity of the phase 
transition, now known as the Griffiths region. However, in classical systems with 
uncorrelated disorder, this Griffiths singularity in the free energy is an essential one 
and thus very weak and probably unobservable in experiment. Disorder correlations 
generically increase the effects of impurities. Therefore, stronger rare region effects 
have been found in classical systems with extended defects and in random quantum 
systems (where the correlations are in imaginary time direction). In the random 
transverse field Ising model [2] , or equivalently, the classical McCoy- Wu model [3] , the 
Griffiths singularity takes a power law form, accompanied by a diverging magnetic 
susceptibility in the Griffiths region. Very recently, it has been found that some phase 
transitions can even be completely destroyed by smearing when the rare regions order 
independently from the bulk system. This happens, e.g., in a classical Ising magnet 
with planar defects [4] and in itinerant quantum ferromagnets [5] . 

In this paper, we investigate the effects of rare regions on nonequilibrium 
phase transitions with quenched spatial disorder. We concentrate on the prominent 
class of absorbing state phase transitions which separate active, fluctuating states 
from inactive, absorbing states where fluctuations cease entirely [6-9]. The generic 
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universality class for absorbing state transitions is directed percolation (DP) [10]. 
According to a conjecture by Janssen and Grassberger [11], all absorbing state 
transitions with a scalar order parameter, short-range interactions, and no extra 
symmetries or conservation laws belong to this class. Examples include the transitions 
in the contact process [12], catalytic reactions [13], interface growth [14], or turbulence 
[15]- 

The effects of uncorrelated spatial disorder, i.e., point-like defects, on the DP 
transition have been studied in some detail in the past. According to the Harris 
criterion [16,17], the DP universality class is unstable against spatial disorder, because 
the (spatial) correlation length exponent v±_ violates the inequality v\_ > 2/d for all 
spatial dimensionalities d < 4. Indeed, in the corresponding field theory, spatial 
disorder leads to runaway flow of the renormalization group (RG) equations [18], 
destroying the DP behavior. Several other studies [19-22] agreed on the instability of 
DP against spatial disorder, but a consistent picture has been slow to evolve. Recently, 
Hooyberghs et al. applied the Hamiltonian formalism [23] to the contact process with 
spatial disorder [24] . Using a version of the Ma-Dasgupta-Hu strong-disorder RG [25] 
these authors showed that the transition (at least for sufficiently strong disorder) is 
controlled by an exotic infinite-randomness fixed point with activated rather than the 
usual power-law scaling. 

Very recently, it has been suggested [26] that extended spatial defects like 
dislocations, disordered layers, or grain boundaries can have an even more dramatic 
effect on noncquilibrium phase transitions in the DP universality class. Rare region 
effects similar to but stronger than the usual Griffiths phenomena [1, 17] actually 
destroy the sharp transition by smearing. This happens because rare strongly coupled 
spatial regions can undergo the transition independently from the bulk system. Based 
on an extremal statistics approach it has been predicted [26] that the spatial density 
distribution in the tail of the smeared transition is very inhomogeneous, with the 
average stationary density and the survival probability depending exponentially on 
the control parameter. 

In the present paper we present results of extensive Monte-Carlo simulations of 
a two-dimensional contact process with linear spatial defects which provide numerical 
evidence for this smearing scenario in a realistic model with short-range couplings. 
The paper is organized as follows. In section we introduce the model and 
briefly summarize the results of the extremal statistics theory for the smeared phase 
transition. In section we present our simulation method and the numerical results 
together with a comparison to the theoretical predictions. We conclude in section 0] 
by discussing the importance of our results and their generality. 

2. Theory 

2.1. Contact process with extended impurities 

The contact process [12] is a prototypical system in the directed percolation 
universality class. It can be interpreted, e.g., as a model for the spreading of a disease. 
The contact process is defined on a d-dimensional hypercubic lattice. Each lattice 
site r can be active (occupied by a particle) or inactive (empty). During the time 
evolution, active sites can infect their neighbors or they can spontaneously become 
inactive. Specifically, particles are created at empty sites at a rate Xn/(2d) where n 
is the number of active nearest neighbor sites and the 'birth rate' A is the control 
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parameter. Particles are annihilated at unit rate. For small birth rate A, annihilation 
dominates, and the absorbing state without any particles is the only steady state 
(inactive phase). For large birth rate A, there is a steady state with finite particle 
density (active phase). The two phases are separated by a nonequilibrium phase 
transition in the DP universality class at A = A°. 

Quenched spatial disorder can be introduced by making the birth rate A a random 
function of the lattice site. Point-like defects are described by spatially uncorrelated 
disorder. We are interested in the case of extended defects which can be described 
by disorder perfectly correlated in d cor dimensions, but uncorrelated in the remaining 
d r = d — d cor dimensions. Here d cor = 1 and 2 corresponds to linear and planar 
defects, respectively. Thus, A is a function of r r which is the projection of the position 
vector r on the uncorrelated directions. For definiteness, we assume that the birthrate 
values A(r r ) are drawn from a binary probability distribution 

P[A(r r )] = (1 -p)<f[A(r r ) - A] + P 6[X(v r ) - cA] (1) 

where p and c are constants between and 1. In other words, there are extended 
impurities of spatial density p where the birth rate A is reduced by a factor c. 

2. 2. Smeared phase transition 

In this subsection, we briefly summarize the arguments leading to the smearing of the 
phase transition and the predictions of the extremal statistics theory [26] to the extent 
necessary for the comparison with the Monte-Carlo results. 

In analogy to the Griffiths phenomena [1,17], there is a small but finite probability 
w for finding a large spatial region of linear size L r (in the uncorrelated directions) 
devoid of impurities. Up to pre-exponential factors, this probability is given by 

w ~ exp(-pLf r ) (2) 

with p = — ln(l — p). These rare regions can be locally in the active phase, even if 
the bulk system is still in the inactive phase. Since the impurities in our system are 
extended, each rare region is infinite in d cor dimensions but finite in the remaining d r 
dimensions. This is a crucial difference to systems with uncorrelated disorder, where 
the rare regions are finite. In our system, each rare region can therefore undergo a 
true phase transition independently of the rest of the system at some X c (L r ) > A°. 
According to finite-size scaling [27], 

X c (L r ) -X° c = AL~+ , (3) 

where <j> is the clean (<i-dimensional) finite-size scaling shift exponent and A is the 
amplitude for the crossover from a d-dimensional bulk system to a 'slab' infinite 
in d cor dimensions but but finite in d r dimensions. If the total dimensionality 
d = d cor + d r < 4, hyperscaling is valid, and <fi = l/v± which we assume from 
now on. 

The resulting global phase transition is very different from a conventional 
continuous phase transition, where a nonzero order parameter develops as a collective 
effect of the entire system, accompanied by a diverging correlation length in 
all directions. In contrast, in our system, the order parameter develops very 
inhomogeneously in space with different parts of the system (i.e., different r r regions) 
ordering independently at different A. Correspondingly, the correlation length in the 
uncorrelated directions remains finite across the transition. This defines a smeared 
transition. 
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In order to determine the global system properties in the vicinity of the smeared 
transition, we combine and @ to obtain the probability for finding a rare region 
which becomes active at A c as 

u;(A c ) ~ exp (-B(A C - A°)- d -^) (4) 

for A c - A° -> 0+. Here, B = pA d ' ^ . 

The total density p (the total number of active sites) at a certain A is obtained by 
summing over all active rare regions, i.e., all regions with A > A c . Since the functional 
dependence on A of the density on any given active island is of power-law type it does 
not enter the leading exponentials but only the pre-exponential factors. Thus, the 
stationary density develops an exponential tail, 

p st (A) ~ exp (-B(A - A°)-<^) , (5) 

for all birth rates A above the clean critical point A^. Analogous arguments can be 
made for the survival probability -P(A) of a single seed site. If the seed site is on an 
active rare region it will survive with a probability that depends on A via a power law. 
If it is not on an active rare region, the seed will die. To exponential accuracy the 
survival probability is thus also given by JSJ. The local spatial density distribution in 
the tail of the smeared transition is very inhomogeneous. On active rare regions, the 
density is of the same order of magnitude as in the clean system. Away from these 
regions it is exponentially small. 

We now turn to the dynamics in the tail of the smeared transition. The long-time 
decay of the density (starting from a state with p = 1) is dominated by the rare regions 
while the bulk contribution decays exponentially. According to finite size scaling [27] , 
the behavior of the correlation time £t of a single rare region of size L r in the vicinity 
of the clean bulk critical point can be modelled by 

&(A, L r ) ~ 4^-^)/^ A - AL- 1 '^ . (6) 

Here A = A — A^ > 0, z is the d-dimensional bulk dynamical critical exponent, and i>± 
and z are the correlation length and dynamical exponents of a d r -dimensional system. 
To exponential accuracy, the time dependence of the total density is given by 

p(X,t) ~ J dL r exp [-pLf - Dt/&(A,L r )) (7) 

where D is a constant. 

Let us first consider the time evolution at the clean critical point A = A°. For 
A = 0, the correlation time © simplifies to £ t ~ L z r . Using the saddle point method 
to evaluate the integral (JJJ, we find the leading long-time decay of the density to be 
given by a stretched exponential, 

lnp(t) ~ -p z /( d r+z) **•/(*•+*) _ (g) 

For A < A°, i.e, in the absorbing phase, the correlation time of the largest islands 
does not diverge but is cut-off by the distance from the clean critical point, £ t ~ A~ ZI/ . 
The large islands with this correlation time dominate the variational integral Q ■ This 
leads to a simple exponential decay with a decay constant r ~ A~ zl/ 

The most interesting case is A > A°, i.e., the tail region of the smeared transition. 
Here, we repeat the saddle point analysis with the full expression © for the correlation 
time. For intermediate times t < t x ~ (A — \0}-( d r+z)v± the decay of the average 
density is still given by the stretched exponential (JBJ. For times larger than the 
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crossover time t x the system realizes that some of the rare regions are in the active 
phase and contribute to a finite steady state density. The approach of the average 
density to this steady state value is characterized by a power-law. 

P (t)-p(oo)~r*. (9) 

The value of tj) cannot be found by our methods since it depends on the neglected 
nonuniversal pre-exponential factors. 

3. Monte-Carlo simulations 

3.1. Simulation method 

We now illustrate the smearing of the phase transition by extensive computer 
simulation results for a 2d contact process with linear defects (d cor = d r = 1). There is 
a number of different ways to actually implement the contact process on the computer 
(all equivalent with respect to the universal behavior). We follow the widely used 
algorithm described, e.g., by Dickman [28]. Runs start at time t = from some 
configuration of occupied and empty sites. Each event consists of randomly selecting 
an occupied site r from a list of all N p occupied sites, selecting a process: creation 
with probability A(r>)/[1 + A(r r )] or annihilation with probability 1/[1 + A(r>)] and, 
for creation, selecting one of the neighboring sites of r. The creation succeeds, if this 
neighbor is empty. The time increment associated with this event is 1/N p . 

Using this algorithm, we have performed simulations for linear system sizes up 
to L = 3000 and impurity concentrations p — 0.2,0.25,0.3,0.35 and 0.4. The relative 
strength of the birth rate on the impurities was c = 0.2 for all simulations. The data 
presented below represent averages of 200 disorder realizations. Because of the large 
system sizes and the high number of realizations, the statistical errors are very low. 
The relative statistical error of the average density ranges from 5p/p ~ 10~ 4 at high 
densities to 8p/p » 10~ 2 for p rj 10~ 3 and to Sp/p » 10 _1 for the lowest shown 
densities p ~ 10~ 6 . Except at these lowest densities, the error is smaller than or 
comparable with the line width in the figures. 

We have chosen the above parameters, viz., a low concentration of impurities 
which have a birth rate much smaller than the bulk, because these conditions are 
favorable for observing the smeared transition in a finite size simulation. If p was too 
large, the exponential drop-offs in eqs. © and I© would be very steep and hard to 
observe over a significant range of A or t, respectively. If c was too close to one, clean 
critical fluctuations would mask the tail of the smeared transition. 

3.2. Time evolution 

In this subsection, we discuss the time evolution of the density starting from a 
completely occupied lattice, p(0) = 1. Figure ^ presents an overview of the behavior 
of a system with impurity concentration p = 0.2, system size L = 3000, and several 
birth rates from A = 1.62 . . . 1.68. The clean critical point is at A^ = 1.6488 [20]. The 
figure shows that the long-time decay of the density in the absorbing phase, A < A°, 
is approximately exponential, in agreement with the expectation discussed after eq. 
iJSJ The decay constant of this exponential increases with decreasing A. In contrast, 
for A > A° the density approaches a nonzero value in the long-time limit. Close to A°, 
the density appears to decay, but slower than exponentially. 
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Figure 1. Overview of the time evolution of the density p for a system with 
L = 3000 and p = 0.2 and several birth rates (A = 1.68, . . . , 1.62 from top to 
bottom) in the vicinity of the clean critical point \® = 1.6488. 




Figure 2. Left: Logarithm of the density at the clean critical point \® 
as a function of t d r/(d r +z) — jO.362 £ or severa ] impurity concentrations (p = 
0.2, ... , 0.4 from top to bottom) and L = 3000. The long-time behavior follows a 
stretched exponential lnp = — _Et ' 362 . Right: Decay constant E of the stretched 
exponential as a function of [— ln(l - p)] z /( d >-+ z ) = [— ln(l -p)] - 638 . 



According to eq. (JSJ, the behavior right at the clean critical point, A = A°, is 
expected to be a stretched exponential rather than a simple exponential decay. To 
shed more light on the time evolution at \ c , the behavior of lnp as a function of 
j-d r /(d r +z) j g p resen ted in the left panel of figure for several impurity concentrations 
p. For our system, d r /(d r + z) = 0.362 with z = 1.76 being the dynamical exponent 
of the clean 2d contact process [29] . The figure shows that the data follow a stretched 
exponential behavior lnp = —Et 362 over more than three orders of magnitude in 
p, in good agreement with eq. (JSJ. (The very slight deviation of the curves from a 
straight line can be attributed to the pre-exponential factors neglected in the extremal 
statistics theory). The right panel of figure |21 shows the decay constant E, i.e., the 
slope of these curves as a function of p = — ln(l — p). In good approximation, the 
values follow the power law E ~ f, z /( d -r+ z ) = ^0.638 p rec ji c ted in JSJ. 

In the tail of the smeared transition, i.e. for A > A° the density has a constant 
nonzero value p sf = p{oo) in the long-time limit. Figure [3] illustrates the approach of 
the density to this value. It shows ln[p(i) — p s t] as a function of ln(t) for several A. 
The long-time behavior is clearly of power-law type, but the exponent depends on A, 
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Figure 3. Double-logarithmic plot of the approach of the density to its nonzero 
stationary value in the tail of the smeared transition for a system with p = 0.2 
and L = 3000 and birth rate A = 1.71, 1.70, 1.69, 1.68 (top to bottom). The long- 
time behavior is of power-law type, (p(t) — p s t) ~ t~ Fits yield exponents of 
approximately 1.00, 1.08, 1.12, and 1.28, respectively. 
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Figure 4. Stationary density p s t as a function of birth rate A for a clean system 
and a system with impurity concentration p = 0.2. System size is L = 1000. 



i.e., it is nonuniversal. These results agree with the corresponding prediction in eq. 
©• 

3.3. Stationary state 

In this subsection we present and analyze the simulation results for the stationary 
state in the tail of the smeared transition, A > A°. Figure 0] shows a comparison of 
the stationary density p st as a function of A between the clean system and a dirty 
system with p = 0.2. The clean system (p = 0) has a sharp phase transition with a 
power-law sing ularity of the density, p st ~ (A - A^)' 3 with f3 w 0.58 in agreement with 
the literature [20]. In contrast, in the dirty system, the density increases much more 
slowly with A after crossing the clean critical point. This suggests either a critical 
point with a very large exponent (3 or exponential behavior. 

Let us now investigate the behavior of the dirty system in the low-density tail 
more closely. In figure [3 we plot lnp st as a function of (A — \ a c )~ drV± for several 
impurity concentrations p, as suggested by eq. (01 . The data in the left panel of 
figure show that the density tail is indeed exponential, following the prediction 
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Figure 5. Left: Logarithm of the stationary density p s t as a function of 
(A — Aj?)~ d rv±_ = (\ _ A^) — 734 for several impurity concentrations p and 
L = 3000. The straight lines are fits to eq. JSJ- Right: Decay constant B as 
a function of — ln(l — p). 



In p s t = —B(X — X° c )~ dTU± over at least two orders of magnitude in p st . (The clean 
2d spatial correlation length exponent is v± = 0.734 [29].) Fits of the data to eq. ((SJ 
are used to determine the decay constants B. The right panel of figure |S] shows these 
decays constants as function of p = — ln(l — p). The dependence is close to linear, as 
predicted below eq. (@J (Slight deviations from the theoretical prediction can again be 
attributed to the pre-exponcntial terms neglected in the extremal statistics theory.) 



4. Discussion and conclusions 



To summarize, we have provided extensive numerical evidence that extended 
impurities destroy the sharp nonequilibrium phase transition in the contact process 
by smearing and lead to a (nonuniversal) exponential dependence of the density and 
other quantities on the control parameter. These results are in agreement with the 
predictions of Ref. [26] which were based on extremal statistics arguments and mean- 
field theory. In this section, we first relate our findings to a power-counting analysis of 
the corresponding field theory. We then compare our smeared phase transition to the 
more conventional Griffiths effects in the contact process with point-like defects [1,17], 
and we discuss general implications for theory and experiment. 

The field-theoretic formulation of the contact process [30, 31] is the so-called 
Reggeon field theory which was originally studied in the context of hadronic 
interactions at ultra-relativistic energies (for a review see, e.g., [32]). In appropriate 
units, its action reads 



S 



d d rdt 



d t 



V 2 + - 



(10) 



where 0(r, t) represents the density while cf>(r, t) is the Martin-Siggia-Rosen response 
field and k denotes the bare distance from the transition. A simple power counting 
analysis for the scale dimension of T at the Gaussian fixed point yields [r] = 4 — d, 
i.e., the upper critical dimension of this field theory is d~£ = 4. Spatially quenched 
disorder (both uncorrelated or correlated) can be taken into account by adding a term 



Sdis — 



dt 



(11) 
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where the outer integral is over the uncorrelated directions r r while the inner spatial 
integrals are over the correlated directions r cor (see [18] for the uncorrelated case 
dcor = 0). Power counting for the scale dimension of 7 gives [7] = 4 — d + d cor . 
Uncorrelated disorder is marginal at d = c£+ = 4, but disorder correlations increase 
the scale dimension of 7 making the disorder term renormalization-group relevant. The 
power-counting analysis thus predicts stronger effects for correlated disorder than for 
point-like disorder, in agreement with our results. Let us emphasize however, that the 
smeared phase transition scenario found in the present manuscript cannot be obtained 
in a perturbative analysis of the Reggeon field theory because the rare regions are non- 
perturbative degrees of freedom. 

Both conventional Griffiths effects and the smearing scenario found in the present 
paper are caused by rare large spatial regions which are locally in the active phase even 
if the bulk system is not. The difference between Griffiths effects and the smearing 
of the transition is the result of disorder correlations. For point-like defects, i.e., 
uncorrelated disorder, the rare regions are of finite size and cannot undergo a true 
phase transition. Instead, they fluctuate slowly which gives rise to Griffiths effects. 
In contrast, if the rare regions are infinite in at least one dimension, a stronger effect 
occurs: each rare region can independently undergo the phase transition and develop 
a nonzero steady state density. This leads to a smearing of the global transition. 

The smearing mechanism found here relies only on the existence of a true phase 
transition on an isolated rare region. It should therefore apply not only to the directed 
percolation universality class, but to an entire family of noncquilibrium universality 
classes for spreading processes and reaction-diffusion systems. Note that while the 
presence or absence of smearing is universal in the sense of critical phenomena (it 
depends on symmetries and dimensionality only), the functional form of the density 
and other observables is not universal, it depends on the details of the disorder 
distribution [4]. 

Smearing phenomena similar to the one found here can also occur at equilibrium 
phase transitions. At quantum phase transitions in itinerant electron systems, even 
point-like impurities can lead to smearing [5] (the necessary disorder correlations are in 
imaginary time direction). In contrast, for the classical Ising (Heisenberg) universality 
class, the impurities have to be at least 2d (3d) for the transition to be smeared which 
makes the phenomenon less likely to be observed [4]. 

In the context of our findings it is worth noting that, despite its ubiquity in 
theory and simulations, clearcut experimental realizations of the directed percolation 
universality class are strangely lacking [33]. To the best of our knowledge, the only 
verification so far has been found in the spatio-temporal intermittency in ferrofluidic 
spikes [34] . We suggest that the disorder-induced smearing found in the present paper 
may explain the striking absence of directed percolation scaling [33] in at least some 
of the experiments. 
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